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Abstract 

In this research it has been proposed to study the effect of filtering processess on range 
images and also to evaluate the performance of two different laser range mappers. 
Median filtering had been utilized to remove noise from the range images. First and 
second order derivatives are then utilized to locate the similarities and dissimilarities 
between the processed and the original images. Range depth information is converted 
into spatial coordinates, and a set of coefficients which describe three dimensional ob- 
jects is generated using the algorithm developed in the second phase of this research. 
Range images of spheres and cylinders are used for experimental purposes. An algo- 
rithm was also developed to compare the performance of two different laser range 
mappers based upon the range depth information of surfaces generated by each of the 
mappers. Further more, an approach based on two-dimensional analytic geometry is 
also proposed which serves as a basis for the recognition of regular three dimensional 
geometric objects. 
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1. Introduction 


The problem of 3-D object recognition has been an interesting research area for 
the past few years with tremendous scope of improvisations in every department of the 
recognition scheme. Unlike the recognition procedures developed for intensity based 
image information, the recent upsurge of several active and passive sensors extracting 
quality range information has lead to the involvement of explicit geometric shapes of 
the objects for the recognition schemes. 

Range images share the same format of the intensity images (i.e. either of these 
images are two dimensional array of numbers), the only difference being that the 
numbers in the range images represent the distances between a sensor focal plane to 
points in space. The laser range finder is the most widely used sensor these days. The 
laser range finder makes use of a laser beam which scans the surfaces in the scene of 
observation from left to right and top to bottom. The distances thus obtained are 
measures of both depth and scanning angle. Until unless a specific algorithm demands 
a special form of these range images, for most of the time it is mainly the depth infor- 
mation which is utlized for the recognition process. 

The range data obtained from a laser radar vision system is chiefly affected with 
two types of problems. The first called the Doppler shift, erupts essentially due to the 
way a laser radar system functions. Recently new radar vision systems have come in 
the market with an inbuilt doppler shift corrector which removes the distortions from 
the range data. The second problem, which is noise in the data picture (mainly salt 
and pepper) is generated on account of the improper wiring circuitry of the whole sys- 
tem. 
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The process by which doppler shift is corrected for our system is discussed in [1]. 
In this report we will be discussing about the median filter which to a large extent 
helps in filtering the noisy range data. 

Median filtering was first suggested by Tukey [5] and since then has been widely 
adopted for two-dimensional image noise smoothing. The most distinguishing property 
of the median filter is that it preserves monotonic step edges, i.e., it does not blur 
sharp edges as most of the linear filters would do. 

Range data from regular objects like spheres, cylinders and cones have been con- 
sidered in this research and the effect of median filtering on each of these has been 
studied. A scheme to evaluate range data obtained from two different laser range 
mappers is also discussed. As the prime objective of this research is to come up with 
a automatic 3-D object classifier, a new approach based upon analytic geometry has 
been proposed for the recognition scheme. 

2. Theoretical Development 

Median Filtering 

Conventionally, a rectangular window of size M x N is used in two dimensional 
median filtering. As in our case, experiments were carried out with square windows of 
mask sizes 3x3 and 5x5. As according to the common belief of the existence of 
salt and pepper at the edges, noise in the range images experimented in this research 
were some what distributed uniformally throughout. Irrespective of the mask size, the 
range information at every pixel in the image is replaced by the median of the the pix- 
els contained in the M x M window centered at that point. Referring to figure 1, keep- 
ing in mind that the dark pixels correspond to the object and the white pixels to the 
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background, specks of white pixels inside the object refers to the salt noise and the 
specks of black pixels in the white background refers to the pepper noise. Figure 3 is 
obtained as a result of a 3 x 3 mask being moved over the entire image. The picture 
looks as sharp as the original image though some of the noise still exists. A 5 x 5 
mask completely removes all the salt and pepper noise, but the image as seen in figure 
4, to some extent has a low contrast, but at the same time has become more smoother 
than the original image. 

Once a range image is filtered using a median filter of different masks, the next 
concern is to study the changes which have been brought about by filtering to the ori- 
ginal data. Evaluating curvatures is one good way of distinguishing similarities and 
dissimilarities among the filtered images and the original range data. 


First and second order derivatives are evaluated along the x- and y-axis to check 
the uniformity of the original and the filtered images. The first order derivative for a 
pixel Aj j centered at i,j is given as: 

3A 1 

0 X ~ 2e [(A ' + i^i — A ij+^ ( A i+ij — A i,j)l> 
and 


0A 1 

~dy ~ 2e^ Ai+1 -i +1 “ A i+i.j) + ( A i.j+i “ A i,j)] 
Similarly the second order derivatives for a pixel centered at Ajj is given as: 


dx 2 e 2 2E ’’ j + W 


and 


0 y 2 £ 2 1 Ai -H 2E i.j + 

e above refers to the spacing between picture cell centers. 
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A sign map whereupon relationship among two neighboring pixels with respect to 
the depth value, is also genertated to make sure that the median filtering does not alter 
the original data to a large extent. 

A second degree general quadric surface as we know is given by the relation, 
F(x,y,z) = ax 2 + by 2 + cz 2 + 2fyz + 2gzx + 2hxy + 2px + 2qy + 2rz + d = 0 

Using the approach formulated by Groshong and Bilbro [1,2] the ten coefficients, 
a, b, c, d, f, g, h, p, q, and r that uniquely describe a quadric surface are determined. 
Coefficients are obtained for each of the filtered images and their relationship with the 
coefficients evaluated for the original range data (one with the noise) are studied for 
each of the surfaces individually. 

Evaluation of the performance of two different laser range mappers. 

In the second phase of our research [1], an approach has been put forward for 
determining the performance of two different laser range mappers using a particular 
test object, i.e., depth maps are obtained for the same object using two different range 
mappers. In this report we have come up with an approach which evaluates the per- 
formance of two different range mappers based upon the depth information obtained 
for two different sizes of the same object, i.e., the test object had the same shape but is 
of different size. The object under consideration is a sphere. 

Theory 

Consider the general equation of the sphere which is in the form of 
(x-x,) 2 + (y- yi ) 2 + (z- Zl ) 2 = 0 (1) 

where xj, y l5 and z } are the coordinates of the center of the sphere. Equation (1) can 
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also be expressed as 

x 2 + y 2 + z 2 + 2fx + 2gf + 2hz + d = 0 (2) 

It is to be noted that the coefficients of x 2 , y 2 , and z 2 are all equal to 1. 

From analytic geometry we know that Xj, y^, and Zj from equation (1) are related 
to the coefficients of x, y, and z in equation (2) with the following relations: 

xj = — 2f 
x 2 = -2g 
and 

x 3 = -2h 

Once the coefficients f, g, and h are evaluated using the algorithm formulated by 
Groshong and Bilbro [2], the center of the sphere, i.e., x, f y h and z x is evaluated using 
the above relationships. It is to be noted that the coefficients f, g, and h and the center 
of the sphere (xj, y x , z x ) evaluated experimentally, certainly do not denote the correct 
coefficients and the center respectively, since a small surface patch of the range data 
has been utilized to determine these coefficients. 

For each set of the sphere range data generated using two different laser range 
mappers, the coordinates of the center of sphere is determined. A least square 
approach as discussed below is next utilized to comment upon the performance of each 
of these laser range mappers. 

Let N be the total number of points (pixels) used to determine the coefficients of 
the sphere generated using laser system 1. 
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Then 


D1 = Z( x r x i) 2 + (yrYi) 2 + (zr z i) 2 

i=0 

where x ; , y i? and Zj are the cartesian coordinates of each of the N depth points, and x lt 
Yi, and z 1 refer to the center of the sphere. 

Now 

j(Dl) 

N 

denotes the mean square error for the system 1. 

A similar approach is carried over for the sphere data generated using sytem 2 
and a mean square error is evaluated. The value of the mean square error determines 
which set of data is more closer to the data generated from a synthetic sphere. 

Object recognition approach based on analytic geometry 

Analytically three dimensional objects are a set of two dimensional curves super- 
imposed upon each other. A sphere for example, is superimposed of circles of varying 
radii. Based upon the 2-D charactersistics of standard curves like circles, parabolas, 
ellipses, and hyperbolas, a unique scheme has been formulated to distinguish standard 
3-D objects like spheres, cylinders, cones and ellipsoids. 

Each object when intercepted with planes in the horizontal and vertical direction 
yields a set of curves which is sufficient enough to recognize each of the objects, and 
at the same time differentiate each from the other. 

Consider the equation of a quadric surface, 

F(x,y,z) = ax 2 + by 2 + cz 2 + 2fyz + 2gzx + 2hxy + 2px + 2qy + 2rz + d = 0 
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If this surface is intercepted with a plane parallel to the yz-axis (which means x is 
a constant), we get a equation of the type 

F(x,y,z) = By 2 + Cz 2 + Fyz+ Qy+ Rz+ D = 0 

which is an equation of a conic. Based upon the discriminat test [4], which says, 

If Ax 2 + Cy 2 + Bxy+ Ex+ Fy+ D = 0 

is a equation of a conic, then, based upon the sign of the discriminant, B 2 -4AC, the 
curves are of three types. 

B 2 - 4AC = 0, 

implies the curve is a parabola. 

B 2 - 4AC < 0, 

implies the curve is an ellipse. 

And finally, 

B 2 - 4AC > 0, 

implies the curve is a hyperbola. 

3. Practical Implementation and Experimental Results 

Two sets of range data namely, the ones generated using system A and system B 
is to be experimented with and the following objectives were to be achieved. Each set 
i.e., A and B are composed of range images of spheres and cylinders respectively. 

1. Study the effect of median filtering of different mask sizes on each of the sets. 

2. Come up with a method which would evaluate the performance of two different 
laser range mappers. 

Making use of the image processing unit in the Image processing and Computer 
Vision lab at ODU, range images of objects like sphere and cylinder were segmented 
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in order to separate the object from the background. 

The resulting image which is referred to as the raw image is then median filtered 
with mask sizes, (a) 3 x 3, (b) 5 x 5, and (c) 7x7. 

Consider figure 1 which is the actual range image of a sphere (belonging to set 
A) with its background. Figure 2 is the image after segmentation. The effect of 
median filtering on figure 2 can be observed in figure 3 (3 x 3 mask), figure 4 (5 x 5 
mask) and figure 5 (7 x 7 mask). 

The curvature sign map which was discussed in the earlier section, is then used to 
study the effect of median filtering on the original image shown in figure 2. Determin- 
ing the first and second derivative with respect to x- and y-axis and comparison of 
each of these maps will determine whether or not the median filtering has altered the 
original range image to any extent. Figures 6(a), 6(b), 6(c), and 6(d) are the first and 
second derivative with respect to x- and y-axis respectively for figure 2. Similarly 
figures 7(a), 7(b), 7(c), 7(d) and figures 8(a), 8(b), 8(c), 8(d) and figures 9(a), 9(b), 
9(c), 9(d) are the first and second derivatives for the figures 3, 4, 5 respectively. 

In all of these figures, the sign "+" is assigned to a particular pixel position if the 
magnitude of the derivative (first or second) of that pixel is greater than the magnitude 
of the derivative (first or second) of the pixel to its right. Similarly the sign is 
assigned to a particular pixel position if the magnitude of the derivative (first or 
second) is lesser than the magnitude of the derivative (first or second) of the pixel to 
its right. In the case when the magnitudes of the derivatives (first and second) of 
either pixels is the same, the sign " " (blank) is assigned. 

Sign maps which were mentioned before are also generated to check the integrity 
of the image data before and after the filtering process. Depending upon the 
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Figure 6(a). First derivative w.r.t x-axis of the original sphere. 
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Figure 6(b). First derivative w.r.t y-axis of the original sphere. 
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Figure 6(c). Second derivative w.r.t x-axis of the original 
cylinder . 
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Figure 6(d). Second derivative w.r.t y-axis of the original sphere. 
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Figure 7(a). First derivative w.r.t x-axis of the sphere filtered 
with a mask size of 3 X 3. 
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Figure 7(b). First derivative w.r.t y-axis of the sphere filtered with 
a mask size of 3 X 3. 
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Figure 7(c). Second derivative w.r.t x-axis of the sphere filtered with 
a mask size of 3 X 3. 
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Figure 7(d). Second derivative w.r.t y-axis of the sphere filtered with 
a mask size of 3 X 3. 
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Figure 8(a). First derivative w.r.t x-axis of the sphere filtered with 
a mask size of 5 X 5. 
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Figure 8(b). First derivative w.r.t y-axis of the sphere filtered 
with a mask size of 5 x 5. 
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Figure 8(c). Second derivative w.r.t x-axis of the sphere filtered with 
a mask size of 5 X 5. 
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Figure 8(d) . Second derivative w.r.t y-axis of the sphere filtered with 
a mask size of 5 X 5. 
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Figure 9(a). First derivative w.r.t x-axis of a sphere filtered 
with a mask size of 7 x 7. 
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Figure 9(b). First derivative w.r.t y-axis of a sphere filtered with 
a mask size of 7 x 7. 
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Figure 9(c). Second derivative w.r.t x-axis of the sphere filtered 
with a mask size of 7 x 7. 
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Figure 9(d). Second derivative w.r.t y-axis of the sphere 
with a mask size of 7 x 7. 
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magnitude of the depth value of a pixel and its adjacent neighbor, a "+” or a sign 
is assigned to the pixel location in the sign map. Figure 10 is the sign map generated 
for the original raw image data of the sphere. Similarly figures 11, 12, and 13 are the 
sign maps for the 3 x 3, the 5 x 5, and the 7 x 7 filtered images of the sphere. A 
careful observation of all these sign maps does suggest that only a small variation has 
been brought about due to the filtering processess. 

Since the main objective of the median filtering is to remove the salt and pepper 
noise in the range images and thus present a noise free range image for the evaluation 
of the object coefficients [1], it is seen from figures 3, 4, and 5 that a fine job has been 
done by all of these filters. However, looking at the curvatures sign maps it is 
observed that, as the mask size of the filter increases, the curvature maps starts looking 
more and more different than the original. The 3 x 3 filtered image being the most 
closest to the original raw image can be utilized for further processing and for 
describing the surface features. 

Once the data files are obtained for each of the images which have been filtered, 
the depth information of each of these files is converted into rectangular coordinate 
system [1]. These cartesian coordinate information is then utilized for determining the 
coefficients which describe each of the objects. 

Listed in table 1 are the coefficients obtained for the original range image , the 3 
x 3 filtered image, the 5 x 5 filtered image and finally the 7 x 7 filtered image of a 
sphere. At a glance none of these coefficient sets for certain describe a real sphere. 
The following procedure is adopted to determine which particular set of coefficients 
best describes the original image data of the object. 
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Figure 10. Sign map generated for the original raw image of the sphere 
taking into consideration the magnitude of the depth value at a 
particular pixel and its neighboring pixel. 
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Figure 11. Sign map generated for the 3x3 filtered image of the sphere 
taking into consideration the magnitude of the depth value at a particular 
pixel and its neighboring pixel. 
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Figure 12. Sign map generated for the 5 x 5 filtered image of the sphere 
taking into consideration the magnitude of the depth value at a particular 
pixel and its neighboring pixel. 
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Figure 13. Sign map generated for the 7 x 7 filtered image of the sphere 
taking into consideration the magnitude of the depth value at a particular 
pixel and its neighboring pixel. 


35 



TABLE 1 


Com 

parison of Coefficients evaluated for the original and the processed images 

Coefficient 

Raw Image 

3x3 filtered image 

5x5 filtered image 

7x7 filtered image 

A, Coeff. of x 2 

0.3026 

0.2211 

-0.4860 

0.4242 

B, Coeff. of y 2 

0.2734 


-0.3291 

0.2178 

C, Coeff. of 2 2 

0.6545 


-0.3338 


E, Coeff. of yz 

0.5310 


0.4834 


F, Coeff. of xz 

0.6357 

-0.4860 

0.7194 


G, Coeff. of xy 

0.3524 

0.2339 

-0.5801 


P, Coeff. of x 

0.30365 

0.19995 


0.3127 

Q, Coeff. of y 

0.4199 

0.4401 

-0.3524 

0.1996 

R, Coeff. of z 

-0.8172 


0.3191 

-0.5858 

D, Constant 

0.2847 

0.3717 

-0.0973 

0.1516 
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A small surface patch of the object is chosen. In the quadratic form 

F(x,y,z) = ax 2 + by 2 + cz 2 + 2fyz + 2gzx + 2hxy + 2px + 2qy + 2rz + d = 0 

the coefficients a, b, c, d, f, g, h, p, q, and r are inserted and for each (x,y,z) of the 
object patch the error is evaluated for each set of coefficients. A plot is thus generated 
in which every point of the surface patch is replaced with the numerals 1, 3, 5, and 7 
signifying that the minimum error was obtained for that particular set of coefficient. 
Numeral 1 refers to the situation when the original set of coefficients fits best, and 
similarly numerals 3, 5, and 7 are used depending whether the 3 x 3 or the 5 x 5 or 
the 7x7 set of coefficients give the least error. Figure 14 is one such plot obtained 
using the coefficients listed in table 1 of the sphere. 

The next objective to achieve is that of evaluating the performance of two 
different laser range mappers. As mentioned before in section 2, the sets A and B 
consist of two different sets of range images abstracted from two different laser range 
mappers. For evaluating the performance, the range information of two different 
spheres obtained from either of these mappers is utilized. Let’s call the range image 
of the sphere using system A as sphere 1. Similarly let’s call the range image of the 
sphere obtained using system B as sphere2. A surface patch of sphere 1 consisting of 
8086 points was selected for experimentation purposes. Similarly the surface patch of 
sphere2 had 726 points. Using the approach discussed in section 2 whereby the mean 
square error is evaluated by trying to a fit a set of data to a real sphere, the mean 
square errors for sphere 1 and sphere2 is obtained. 

Mean square errors are obtained for the raw image, and the 3 x 3 image for 
sphere 1 and sphere2. The mean square error for the sphere 1 belonging to set A was 
found to be 0.010191 units and 0.009921 units (raw image and 3x3 filtered image 
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Figure 14. Best fit plot for the sphere belonging to set A. Numerals "1, 3, and 5" denote 
the original sphere, 3x3 filtered image, and 5x5 filtered image respectively. 
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respectively). The mean square errors for sphere2 belonging to set B was found to be 
.019095 units and 0.018686 units (raw and 3x3 filtered images respectively). 

The curvature maps for sphere and cylinder belonging to the sets A and B are 
shown in appendix A. Appendix B lists out the ten coefficients obtained for all the 
different images of sets A and B. Files with extension *.cod serve as the input for the 
program evaluating the coefficients, and the files with extension *.coe consists of the 
output data, which are the needed necessary coefficients. Appendix C consists of a 
detailed listing of all the programs utilized. 

4. CONCLUSIONS 

In this research, range images of objects obtained using laser range mappers are 
utilized to recognize three dimensional regular objects. Due to inherent problems in 
the laser range mappers, the depth information obtained by itself cannot be utilized to 
make a accurate description of the object. The approach involving the evaluations of 
the ten coefficients which best describe an object is utilized on filtered images of the 
original objects. Inspite of using noise free images, it is seen that the coefficients 
obtained for each object does not infer the shape of any of the objects. 

A new approach which involves 2-D analytical geometry has been discussed 
briefly which appears very promising for the recognition of 3-D objects. The 
coefficients obtained earlier do come in handy while using a discriminant test for 
describing each of the objects with 2-D curves. In the future research the above new 
theory formulated will be utilized for making a accurate description of each of the reg- 
ular 3-D objects. 

Calculations evaluating the performance of the two different laser range mappers 
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quite distinctly showed that laser range mapper A performs better than laser range 
mapper B. 
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APPENDIX A 

Curvature sign maps of the following range images is included in this appendix. 

1. Original cylinder image belonging to set A. 

2. 3x3 filtered image of the cylinder belonging to set A. 

3. 5x5 filtered image of the cylinder belonging to set A. 

4. 7x7 filtered image of the cylinder belonging to set A. 

5. Original sphere image belonging to set B. 

6. 3x3 filtered image of the sphere belonging to set B. 

7. 5x5 filtered image of the sphere belonging to set B. 

8. Original cylinder image belonging to set B. 

9. 3x3 filtered image of the cylinder belonging to set B. 

10. 5x5 filtered image of the cylinder belonging to set B. 

For each of the above images the curvature sign maps consists of the first and 
second derivatives with respect to the x- and y-axis. Sets A and B signify to the 
images mapped by two different laser range mappers. 
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Images belonging to set A 
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++++++++++-> f - H -H-H-f fr-H - 
-H-+-H-++++-H-+ lilt )■+ + +++ 
- +++++++++++ + +++++ - +++ 
+ - ++++++++++++++ - - +++ - 

++ +++++ I I t <■■« + 

+++++++ +++++++ - - + 

- +++++++++++++++++++++ 
+ - - ++++++++++++ - +++++ - 

+++ ++++ +++++ - + 

+-H-++++ ++++++++ - ++ 


-H- - +++++++++++++ - ++++ - 
+++++++++++++++++++-H-+ 
- - +++++++++++++++ - - ++- 
++ - - ++++++++++ +++ - + 

- ++-H-++++++++++++++ - + - 

+ - - +++++++++++++ +++ 

-H — h+H — I — t~H — I — ! — - - - 

- -++ -H- + - + 


First derivative w.r.t y-axis of the cylinder filtered with 
a mask size of 3 X 3. 
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+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 


Second derivative w.r.t x-axis of the cylinder 
filtered with a mask size of 3 X 3 . 
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+++ I I I M H - 

t I- +++++++++4 t I t - 

++- 

I II II II I +++++-H-H-++- - + 
- +++++++++++++++++ - +++ 
+ - I I )+++++++++++- -+++- 


- + +++++ ++++++ 

+ -+++--- 

- ++++++++-H-+++-H- +++ 

I l l + ++++++++ - +++++ - 

+++ ++++ +++++ - + 


M"4"f+-M-+++++++ - I M-+ 

+++++++++++++ ++ - 

+++ +++- - + 

- - +++ - +++++++++++ - -++- 

+ ++++++++++ + - + 

-+++ +++++- + - 

+ ++ +++ 

+++ - ++++++++++ - - +++ - + - 


Second derivative w.r.t y-axis of the cylinder 
filtered with a mask size of 3 X 3. 
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++++++++++++-> i H H H 4 - ++ 


+ - -+++-H-+-HH--M-++- 

+++ H H- + ++++- + 

+++++++ f f i - HH44 I H M H 
- ++- H -+++4H-+-HH-+++ - - +++ 

+ - - +++++++++++ ++++- 

+++ ++++ ++++++ - + 


- +++++++++++ + + +++ - -++- 
+ +++++++++++ - - +++++ 


- -+4HH-+++++^H-++++- -++- 
+-f - - ++++++-H-++ +++ - + 


- ^+-HH-+++++++++++ - - + - 

H HH-++++++++-H ++ - + 

++++++++ I 1+-H HH f 4 t + - + - 

- -+++ ++ + - + 


First derivative w.r.t y-axis of the cylinder 
filtered with a mask size of 5 X 5. 
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- -H-++++++++++ - - 



First derivative w.r.t x-axis of the cylinder 
filtered with a mask size of 5 X 5 . 
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+ 


+ 


+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
-+ 
+ 
■+ 
*+ 
■+ 
■ + 
■ + 
■ + 


Second derivative w.r.t x-axis of the cylinder filtered with 
a mask size of 5 X 5. 
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+++ ++++++++ 

- ++++++++++++++++++ - ++ 


++++++++++++++ - ++ 
- -H-H-++++++++-H-H-+ - +++ 
+ - - ++-H-++++-H-+++ - - +++ - 
+++++ + - - 


+++++++ 

+++++++++++ ++++ - 

++++ + - - 

++++++H ++++++++ - ++ 

- ++++++ +++++++++ - - ++ - 

+++++++++++ ++ - 

++++ ++++-- + 

— H+++++ I (■ +++ - ~ ++ - 

+ ++++++++++ ++ - + 

-+++ ++++-- + - 

- - -++ ++ + - + 

+++- -+++++++++- -+++-+- 
- +++++++++++++++++++++ 


Second derivative w.r.t y-axis of the cylinder filtered with 
a mask size of 5 X 5. 
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APPENDIX B 

This appendix consists of the ten coefficients generated for the original and pro- 
cessed range images of a sphere and cylinder mapped using two different laser range 
mappers. Files with extension *.cod consists of range data converted into cartesian 
coordinates, and the files with extension *.coe consists of the coefficients generated for 
each of the images. 
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The input file was "sprawl.cod 
The output file is "sprawl. coe 
The coeff of x-squared is 0.3026157 
The coeff of y-squared is 0.2734349 
The coeff of z-squared is 0.6545654 


The coeff of 

yz 

is -0.5310194 

The coeff of 

zx 

is -0.6357662 

The coeff of 

xy 

is 0.3524517 

The coeff of 

X 

is 0.3036514 

The coeff of 

y 

is 0.4199182 

The coeff of 

z 

is -0.8172019 

The constant 

d 

is 0.2847408 

Coefficients of the 

original sphere image belonging to group A. 
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The input file was "spraw31.cod 
The output file is "spraw31.coe 
The coeff of x-squared is 0.2211579 
The coeff of y-squared is 0.2802473 
The coeff of z-squared is 0.7747064 
The coeff of yz is -0.5038247 

The coeff of zx is -0.4860164 

The coeff of xy is 0.2339016 

The coeff of x is 0.1995363 

The coeff of y is 0.4401489 

The coeff of z is -1.016356 

The constant d is 0.3717703 


Coefficients of the 3 x 3 filtered image of the sphere belonging to group A. 
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The input file was "spraw51.COD 
The output file is "spraw51.COE " 
The coeff of x-squared is -0.4860452 
The coeff of y-squared is -0.3291118 
The coeff of z-squared is -0.3338964 
The coeff of yz is 0.4834592 

The coeff of zx is 0.7194569 

The coeff of xy is -0.5801437 

The coeff of x is -0.3159497 

The coeff of y is -0.3524498 

The coeff of z is 0.3191445 

The constant d is -9.7348504E-02 


Coefficients of the 5x5 filtered image of the sphere belonging to group A. 
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The input file was "sprawmel.cod 
The output file is "sprawmel.coe 
The coeff of x-squared is 0.4242373 
The coeff of y-squared is 0.2178874 
The coeff of z-squared is 0.5845248 
The coeff of yz is -0.3417171 

The coeff of zx is -0.7452961 

The coeff of xy is 0.4353395 

The coeff of x is 0.3127908 

The coeff of y is 0.1996729 

The coeff of z is -0.5858592 

The constant d is 0.1516084 


Coefficients of the 7 x 7 filtered image of the sphere belonging to group A. 
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The input file was "cyrawl.cod 
The output file is "cyrawl.coe " 
The coeff of x-squared is 0.1555596 
The coeff of y-squared is 0.2353804 
The coeff of z-squared is 0.8288453 


The coeff of 

yz 

is 

-0.6818960 

The coeff of 

zx 

is 

3.70348 17E-02 

The coeff of 

xy 

is 

2.1725880E-02 

The coeff of 

X 

is 

-0.2105054 

The coeff of 

y 

is 

0.5823037 

The coeff of 

z 

is 

-1.317142 

The constant 

d 

is 

0.5681907 


Coefficients of the original cylinder belonging to group A. 
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The 
The 
The 
The 
The 

The coeff of 
The coeff of 
The coeff of 
The coeff of 
The coeff of 
The coeff of 
The constant d 


is -0.5485628 
is 0.5481051 
is -0.2466192 
is -0.7515414 
is 0.5662742 
is -1.360964 
is 0.6880789 


input file was "cyraw31.cod 
output file is "cyraw31.coe " 
coeff of x-squared is 0.2676638 
coeff of y-squared is 0.1930158 
coeff of z-squared is 0.7483451 
yz 
zx 
xy 
x 

y 

z 


Coefficients of the 3x3 filtered image of the cylinder belonging to group A. 
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The input file was "cyraw51.cod 
The output file is "cyraw51.coe 
The coeff of x-squared is 5.4338872E-02 
The coeff of y-squared is 9.9206299E-02 
The coeff of z-squared is 0.2060992 
The coeff of yz is -0.1109364 

The coeff of zx is 1.265334 

The coeff of xy is -0.5254330 

The coeff of x is -1.185869 

The coeff of y is 0.3039300 

The coeff of z is -0.7311586 

The constant d is 0.5089003 


Coefficients of the 5x5 filtered image of the cylinder belonging to group A. 
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The input file was "cyrawmel.cod 
The output file is "cyrawmel.coe 
The coeff of x-squared is 0.1532317 
The coeff of y-squared is -9.9520542E-02 
The coeff of z-squared is -0.4889523 


The coeff of 

yz 

is 

0.4767834 

The coeff of 

zx 

is 

1.008621 

The coeff of 

xy 

is 

-0.4587431 

The coeff of 

X 

is 

-1.006533 

The coeff of 

y 

is 

-0.2328676 

The coeff of 

z 

is 

0.4734453 

The constant 

nn* • 

d 

is 

-1.3768099E-02 


Coefficients of the 7 x 7 filtered image of the cylinder belonging to group A. 
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The input file was "R3SPHERE.COD 
The output file is "R3SPHERE.COE 
The coeff of x-squared is 0.1027336 
The coeff of y-squared is 3.8939383E-02 
The coeff of z-squared is 0.5696317 
The coeff of yz is 0.6472183 

The coeff of zx is -0.9516000 


The coeff of xy 

The coeff of x 

The coeff of y 

The coeff of z 

The constant d 


is -4.9645 115E-02 
is 0.8889613 
is -0.6169493 
is -1.387174 
is 0.7926015 


it 


Coefficients of the original sphere belonging to group B. 



The input file was "R3SPHER3.COD 
The output file is "R3SPHER3.COE 
The coeff of x-squared is -4.9067583E-02 
The coeff of y-squared is 0.4412566 
The coeff of z-squared is 0.6636547 
The coeff of yz is -1.2313786E-02 

The coeff of zx is -0.5175490 

The coeff of xy is -0.6759338 

The coeff of x is 0.6078625 

The coeff of y is -0.3856263 

The coeff of z is -1.276605 

The constant d is 0.6796699 


Coefficients of the 3x3 filtered image of the sphere belonging to group B. 
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The input file was "R3SPHER5.COD 
The output file is "R3SPHER5.COE " 
The coeff of x-squared is -5.7173960E-02 
The coeff of y-squared is -0.1170360 
The coeff of z-squared is 0.5475225 
The coeff of yz is 0.5604561 

The coeff of zx is 1.006907 

The coeff of xy is 0.1962991 

The coeff of x is -1.006644 

The coeff of y is -0.3863406 

The coeff of z is -0.9040802 

The constant d is 0.3384019 


Coefficients of the 5 x 5 filtered image of the sphere belonging to group B. 
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The input file was "R6CYLIN.C0D 
The output file is "R6CYLIN.COE 
The coeff of x-squared is 0.9754460 
The coeff of y-squared is 2.5 132844 E-02 
The coeff of z-squared is 3.5924029E-02 


The coeff of 

yz 

is 

-6.8559073E-02 

The coeff of 

zx 

is 

3.1578626E-02 

The coeff of 

xy 

is 

0.2957501 

The coeff of 

X 

is 

0.2924450 

The coeff of 

y 

is 

0.1052131 

The coeff of 

z 

is 

-1.9418295E-02 

The constant 

d 

is 

1.5252778E-02 

Coefficients of 

the 

original cylinder belonging to group B. 
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The input file was "R6CYLIN3.COD 
The output file is "R6CYLIN3.COE 
The coeff of x-squared is -4.7388867E-02 
The coeff of y-squared is -0.3104874 
The coeff of z-squared is -0.3682815 
The coeff of yz is 1.192302 

The coeff of zx is 0.1264399 

The coeff of xy is -0.3063811 

The coeff of x is -2.7492255E-02 

The coeff of y is -0.9607195 

The coeff of z is 0.3220469 

The constant d is 4.0601 194E-03 


Coefficients of the 3 x 3 cylinder image belonging to group B. 
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The input file was "R6CYLIN5.COD 
The output file is "R6CYLIN5.COE 
The coeff of x-squared is 1.76 1973 IE-02 
The coeff of y-squared is 0.7016529 
The coeff of z-squared is -0.2045088 
The coeff of yz is -0.3910733 

The coeff of zx is -0.7922655 

The coeff of xy is -0.3879120 

The coeff of x is 0.8651381 

The coeff of y is -0.1430389 

The coeff of z is 0.2737453 

The constant d is 1.2079749E-02 


Coefficients of the 5x5 filtered image of the cylinder belonging to group B. 
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APPENDIX C 


This appendix consists of the listings of the following programs. 

1. Program which performs the 3 x 3 and 5x5 median filtering. 

2. Program that evaluates the first and second derivative w.r.t to x- and y-axis of the 
data files and then transforms it into a sign map. 

3. Program that generates the sign map for each of the range images based upon the 
magnitude of range value of neighboring pixels. Sign maps for the cylinder of 
set A and the sphere and cylinder of set B are included at the end of the listing. 

4. Program that generates a numeral map based upon the evaluation of the least 
square errors from the generated coefficients. 

5. Program that generates the ten coefficients which describes each of the range 
images. 
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C**** PROGRAM MEDIAN FILTERING 


C**** THIS PROGRAM PERFORMS THE MEDIAN FILTERING ON THE 
C**** ORIGINAL RANGE IMAGE FILES. BY CHANGING THE 
C’*** PARAMETER "M". A 3x3 OR A 5x5 MASK SIZE CAN BE UTILIZED 
C**** F or FILTERING. 


PARAMETER (N=512) 

INTEGER*2 A(N,N),MED(N,N) 

CHARACTER *12 INFILE,OUTFILE 

C 

C MAIN PROGRAM 

C 

WRITE(*,123) 

123 FORMAT (5X,’INPUT FILE NAME : INFILE’) 

READ(*,*)INFTLE 
WRITE(*,223) 

223 FORMAT(5X,’OUTPUT FILENAME : OUTFTLE’) 

READ(*,*)OUTFILE 


OPEN (UNIT = l.FILE = INFILE, RECL = 2048,STATUS = ’OLD’) 
READ (l,9)((A(yy = 1,N),I = 1,N) 

9 FORMAT(512I4) 

M=3 

CLOSE(l, DISPOSE = ’SAVE’) 

CALL MEDFLT(A,MED,N,M) 

OPEN (UNIT = 2,FILE = OUTFILE.RECL = 2048.STATUS = ’NEW’) 

WRITE (2,ll)((MED(y)J = l,N),I = l,N) 

11 FORMAT (5 1214) 

CLOSE(2,DISPOSE = ’SAVE’) 

STOP 

END 

cc 

CC SUBROUTINE MEDIAN FILTER 

CC 

SUBROUTINE MEDFLT (A,MED,N,M) 

INTEGER*2 A(N,N),MED(N,N),SORT(50) 

LOGICAL NEXCHAN 
C 
C 
C 

MM=M ** 2 
X=(M + l)/2 
Y=X-1 

Ml=(MM + l)/2 
DO 7 I=X,(N-Y) 

DO 9 J=X,(N-Y) 

K1=0 

DO 11 K=(I-Y),(I + Y) 
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'J VO 


DO 13 L=(J-Y),(J + Y) 

K1 = K1 + 1 
SORT(Kl) = A(K,L) 

13 CONTINUE 

11 CONTINUE 

DO 15 U = 1,(MM-1) 

DO 17 K1 = 1,(MM-I1) 

IF (SORT(Kl).GT.SORT(Kl + 1)) THEN 
TEMP = SORT(Kl) 

SORT(Kl) = SORT(Kl + 1) 

SORT(Kl + 1)=TEMP 
END IF 
17 CONTINUE 

15 CONTINUE 

MED(U)=SORT(Ml) 

CONTINUE 
CONTINUE 
DO 19 I = 1,Y 
DO 21 J = 1,N 
MED(y)=A(U) 

MED(N + 1-U) = A(N + 1-U) 

MED(J,N + 1-1) = A(J,N + 1-1) 

MED(J,I) = A(J,I) 

21 CONTINUE 

19 CONTINUE 

RETURN 
END 
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C***** PROGRAM DERIVATIVES 


THIS PROGRAM DETERMINES THE DERIVATIVES 
C***** ALONG THE X-AXIS AND THE Y-AXIS. A GROUP OF FILES CAN BE 

COMAPARED TO SEE WHETHER A PARTICULAR LOCATION HAS THE SAME 
C*.**. CURVATURE OR NOT. 


20 


INTEGER *2 
REAL 
REAL 
REAL 
REAL 
INTEGER*2 
CHARACTER *12 
CHARACTER *2 


IU1,T1,P1,K,L,U 

DX1,DX2,DX3,D Y 1,DY2,DY3 

DX11,DX22,DX33,DY11,DY22,DY33 

D(70350),E(70,350)^(1000,3)A^(60,50) 

D1(70,350),E1(70,350) 

STREC.ENDREC 
INFILE1,INFILE2,INFILE3, POINT 
GRAPH1(70,100),GRAPH2(70,100),GRAPH3(70 ) 100) 


CHARACTER *2 GRAPH4(70,100) 

WRITE(*,20) 

FORMAT(5X, ’INPUT FILE NAME : INFILET) 
READ(*,*)INFILE1 

OPEN(UNTT « 1, FILE = INFTLE1, STATUS = ’UNKNOWN*) 


WRITER, 25) 

25 FORMAT (5X,’INPUT TOTAL # OF PTS : Nl’) 

READ(*,*)N1 
DO 100 I = 1,N1 
READ(l,*)(A(I,jy=l > 3) 

100 CONTINUE 

DO 811 K=l,51 
DO 815 L=l,19 
AA(K,L)=A(L+(19*(K-1))3) 

815 CONTINUE 
811 CONTINUE 


300 FORMAT(512I4) 

C** TO FIND THE DERIVATIVE ALONG X-AXIS 
Cllll WRITE(*,908) 

C908 FORMAT (’INPUT THE STARTING RECORD NUMBER: STREC) 

C READ(*,*)STREC 
C WRITER, 9008) 

C9008 FORMAT (’INPUT THE ENDING RECORD NUMBER: ENDREC’) 

C READ(\*)ENDREC 

OPEN(UNIT = 2,FILE = ’FILE1 X\ STATUS = ’UNKNOWN*) 

OPEN(UNIT = 3, FILE = ’FILEl.Y’.STATUS = ’UNKNOWN’) 

OPEN(UNIT = 4,FILE = ’FILE1JOC’, STATUS = ’UNKNOWN’) 
OPEN(UNIT = 8 JILE = ’FTLE1.YY’, STATUS = ’UNKNOWN’) 

11178 DO 1104 11=1,51 
DO 1204 Jl=l,19 

D(iy 1) = O^^CAACiyi + l)-AA(IUl)) + (AA(I1 + Ml + 1)-AA(I1 + 1J1))) 
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D1(I1,J1) = (AA(I1JM)-2*(AA(IU1» + AA(IM1 + 1)) 


E1(I1J1) = (AA(I1 + U1)-2*(AA(I1J1)) + AA(Il-Ul)) 

E(IU1) = 04*((AA(I1,J1 + 1)-AA(IU1 + 1)) + (AA(Il r Fl)-AA(Il + VI))) 
1204 CONTINUE 
1104 CONTINUE 
1965 DO 11104 11 = 1,51 

WRITE(2,*)(D(IUiy 1 = 1,19) 

WRITE(3,*)(E(IUl) r Il = 1,19) 

WRITE(4,*)(Dl(IUiyi= 1,19) 

WRITE(8,*)(El(Il r Jl) r Il= 1,19) 

11104 CONTINUE 
CLOSE(2) 

CLOSE(3) 

CLOSE(4) 

CLOSE(8) 

OPEN(UNIT = 2, FILE = ’FILE1.X’, STATUS = ’UNKNOWN’) 
OPEN(UNTT = 3, FILE = ’FILE1.Y*, STATUS = ’UNKNOWN*) 
OPEN(UNIT = 4, FILE = ’FILE1.XX’, STATUS = ’UNKNOWN*) 
OPEN(UNIT = 5.FILE = ’FILE1.YY’, STATUS = ’UNKNOWN’) 

DO 324 11 = 1,51,1 
READ(2,*)(D(Il r Il) r Jl = 1,19) 

324 CONTINUE 
DO 325 11 = 1,51,1 
DO 326 Jl= 1,19 

IF (D(iyi).GT.D(Il r II + 1))THEN 
GRAPHl(iyi)= ’-’ 

ENDIF 

IF (D(iyi).LT.D(iyi+ 1))THEN 
GRAPH 1(114 1) = ’ + ’ 

ENDIF 

IF (D(iyi).EQ.D(iyi+ 1))THEN 
GRAPHl(iyi) = ’ ’ 

ENDIF 

326 CONTINUE 

325 CONTINUE 
DO 328 11 = 141,1 
READ(3,*)(D1(IU1),J 1 = 1,19) 

328 CONTINUE 

DO 329 11 = 141,1 
DO 330 J1 = 1 19 

IF (D1(IU1).GT,D1(I14I+1))THEN 
GRAPH2(Il r Il) = ’-’ 

ENDIF 

IF (D1(IU1).LT.D1(I14I+ 1))THEN 
GRAPH2(I141)= ’+’ 

ENDIF 

IF (Dl(iyi).EQ.Dl(Il r n+l))THEN 
GRAPH2(iyi)= ’ ’ 
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ENDIF 

330 CONTINUE 
329 CONTINUE 

DO 33211=1,51,1 
READ(4,*)(E(Il r Il) r Jl = 1,19) 

332 CONTINUE 
DO 333 11 = 1,51,1 
DO 334 Jl = 1,19 

IF (E(U r Il).GT.E(Il r II+ 1))THEN 
GRAPH3<Il r ri)= *-* 

ENDIF 

IF (E(I1J1).LT.E(IUI+ 1))THEN 
GRAPH3(Il r Fl)= ’ + ’ 

ENDIF 

IF (E(Il r Jl).EQ.E(Il r II + 1))THEN 
GRAPH3(IU1)= ’ ’ 

ENDIF 

334 CONTINUE 

333 CONTINUE 

DO 336 11 = 1,51,1 
READ(5,*)(El(Il r Il) r Il = 1,19) 

336 CONTINUE 
DO 337 11 = 1,51,1 
DO 338 Jl = 1,19 

IF (E1(I1,J1).GT.E1(I1,JI + 1))THEN 
GRAPH4(IU1)= 

ENDIF 

IF (E1(I1,J1).LT.E1(I1,JI + 1))THEN 
GRAPH4(Il r Il)= ’+’ 

ENDIF 

IF (El(Il r n).EQ.El(Il r JI + l))THEN 
GRAPH4{Il r Il) = ’ ’ 

ENDIF 

338 CONTINUE 

337 CONTINUE 
1324 CONTINUE 

OPEN(UNIT = 13,FTLE = ’GRAPH X\ STATUS = ’UNKNOWN’) 
OPEN(UNIT = 14, FILE = ’GRAPH.Y’, STATUS = ’UNKNOWN’) 
OPEN(UNIT = 15, FILE = ’GRAPH .XX’, STATUS = ’UNKNOWN’) 
OPEN(UNIT = 16, FILE = ’GRAPH. YY’, STATUS = ’UNKNOWN’) 
DO 21104 11=141,1 

WRITE(13,1234)(GRAPHl(Il r Il) r Il = 1,19) 
WRITE(14,1234)(GRAPH2(Il r Il) r Il = 1,19) 
WRITE(15,1234)(GRAPH3(Il r Jl) r Jl = 1,19) 
WRITE(16,1234)(GRAPH4(Il r Il) Jl = 1,19) 

21104 CONTINUE 
1234 FORMAT (30X,2QA1) 

C WRITE(*,21) 

C GOTO 64 

END 
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C****PROGRAM RANGE SIGN MAP 


C***** THIS PROGRAM GENERATES A SIGN MAP FOR DATA FILES BY TAKING 
INTO CONSIDERATION THE ABSOLUTE DIFFERENCE IN RANGE VALUE 
C***** OF NEIGHBORING PIXELS. 


INTEGER *2 A(0:5U,0:512),D(100,100) 

INTEGER*2 U,J1,T1,P1,ZZ^X 

CHARACTER * 12 INFILE1,INFILE2,INFILE3, POINT 

CHARACTER *2 GRAPH1(100,100) 

WRITE(*,20) 

20 FORMAT(5X, ’INPUT FILE NAME : INFILE1’) 

READ(*,*)INFILE1 

OPEN(UNIT = 1, FILE = INFILE 1, STATUS = ’UNKNOWN’, RECL=2048) 
DO 100 1 = 1,511 
READ(l,300)(A(Uy = 1,512) 

100 CONTINUE 
300 FORMAT (5 1214) 

ZZ=1 
C XX = 1 

DO 43 1 = 165,215 
XX= 1 

DO 53 J = 260,278 
D(ZZ^X)=A(U) 

C ZZ=ZZ+1 

XX=XX+1 
53 CONTINUE 

C XX=1 

ZZ=ZZ+1 
C XX= 1 

43 CONTINUE 

WRITE(V)XX,ZZ 


C**** test file used for this program is that of the cylinder 

C**** BELONGING TO SET A. 

OPEN(UNIT = 2.FTLE = ’RANGEVAL.DAr .STATUS = ’UNKNOWN) 
OPEN(UNIT = 3, FILE = ’RANGEDIFF.DAT, STATUS = UNKNOWN) 
c OPEN(UNIT = 4, FILE = ’FILE1XX’, STATUS = ’UNKNOWN’) 


DO 325 I=1,ZZ-1 
DO 326 J = l.XX-1 
IF (D(U).GT.D(U + 1))THEN 
GRAPH1(U)= ’ + ’ 

ENDIF 

IF (D(U).LT.D(U + 1))THEN 
GRAPH1(U)= 
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ENDIF 

IF (D(U).EQ.D(y+l))THEN 
GRAPH l(y)= ’ * 

ENDIF 

326 CONTINUE 
325 CONTINUE 

DO 21104 1 = 1,ZZ-1 

WRITE(3,1234)(GRAPHl(I r J)rf = 1,XX-1) 
WRrrE(23000)(D(U),J = 1,XX-1) 

21104 CONTINUE 
1234 FORMAT (35X,20A 1) 

3000 FORMAT(I4) 

STOP 

END 


107 



c****program best fit coefficients 


C***« THIS PROGRAM MAKES A PLOT USING THE COEFFICIENTS GENERATED 
C*«. FROM THE PROGRAM "SURFACE.FOR". AT EACH PIXEL OF A TEST 
C***» SURFACE PATCH, THE ERROR IS DETERMINED USING THE GENERATED 
C**** COEFFICIENTS OF THE ORIGINAL RANGE DATA, THE 3X3 RANGE IMAGE, 
C**** THE 5X5 RANGE IMAGE, AND THE 7X7 RANGE IMAGE. WHICHEVER 
C***’ GIVES THE MINIMUM ERROR REPLACES THE PIXEL WITH THE NUMERAL 
C**** 1, 3, 5, 7 WHEREEVER APPLICABLE. 


REAL A(5000,3),B(5000,3),C(50003),D(5000),H(5000,3) 

REAL E(5000),F(5000),P(5000) 

INTEGER G(5000),PLOT (100, 100) 


C**** TEST FILE IN THE PROGRAM ARE THE RANGE IMAGES OF THE 
C**** CYLINDER BELINGING TO GROUP A. 

OPEN(UNTT = 1,FILE = ’CYRAW1.PLT, STATUS = ’UNKNOWN’) 
OPEN(UNIT = 2, FILE = ’CYRA WME1.PLT\ STATUS = ’UNKNOWN’) 
OPEN(UNIT = 3, FILE = ’CYRAW51.PLT\STATUS = ’UNKNOWN’) 
OPEN(UNIT = 4, FILE = ’CYRAW31.PLT’, STATUS = ’UNKNOWN’) 
OPEN(UNIT = 8, FILE = ’CYLINDE2.PLT’,STATUS = ’UNKNOWN’) 

DO 10 1 = 1,969 
READ(l,*)(A(Uy=U) 

10 CONTINUE 

DO 40 1=1,969 
C DO 50 J = 1,3 

D(I) = (0.15555 *A(1, 1) * A(1, 1)) + (,23538*A(I,2)*A(I,2)) + 

+ (0.8288*A(I,3)*A(I,3))-(0.6818*A(I,2)*A(I3)) + 

+ (0.03703*A(I,1)*A(I,3)) + (0.021725*A(I,1)*A(I,2))- 

+ (0.2105 *A(I,1)) + (038230* A(1, 2))- 

+ (1.317142*A(I3)) + (0.568190) 

40 CONTINUE 
DO 20 1 = 1,969 
READ(2,*)(B(Uy=l,3) 

20 CONTINUE 

DO 50 1 = 1,969 
C DO 50 J = 13 

E(I) = (0.15323*B(I,1)*B(I,1))-(.09952*B(I,2)*B(I,2))- 
+ (0.48895*B(I3)*B(I3)) + (0.47678*B(I,2)*B(I,3)) + 

+ (1.00862*B(I,1)*B(I3))-(0.4587431*B(I,1)*B(I,2))- 

+ (1.006533*B(I,1))-(0.23286*B(I,2)) + 

+ (0.473445*B(I3))-(0.013768) 

50 CONTINUE 

DO 30 1=1,969 
READ(3,*)(C(U),J=13) 

30 CONTINUE 

DO 60 1=1,969 
C DO 50 J = 13 
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F(I) = (0.054338*C(I,1)*C(I,1)) + (.099206*C(I,2)*C(I,2)) + 
+ (0.2060992*C(I,3)*C(I,3))-(0.110936*C(I,2)*C(I,3)) + 

+ (1.265334*C(I,1)*C(U))-(0.525433*C(I,1)*C(I,2))- 

+ (1.18586*C(I,1)) + (0303930*C(I,2))- 

+ (0.7311586*0(1,3)) + (0.5089003) 

60 CONTINUE 

DO 301 1 = 1,969 
READ(4,*)(H(W=U) 

301 CONTINUE 

DO 602 I = 1,969 
C DO 50 J = 1,3 

P(I) = (0.26766*H(I,1)*H(I,1)) + (.193015*H(I,2)*H(I,2» + 
+ (0.7483451*H(I,3)*H(I,3))-(0.548105*H(I,2)*H(I,3)) + 

+ (0.548105*H(I,1)*H(I,3))-(0.246619*H(I,1)*H(I,2))- 

+ (0.751541*H(I,1)) + (0.5662742*H(I,2))- 

+ (1360964*H(I,3)) + (0.6880789) 

602 CONTINUE 

DO 90 I = 1,969 

IF((D(I).LT.E(I))AND.(D(I).LT.F(I))AND. 

+ (D(I).LT.P(I)))THEN 
G(I) = 1 
ENDIF 
C ENDIF 

IF((E(I).LT.D(I))AND.(E(I).LT.F(I))AND. 

+ (E(I).LT.P(I)))THEN 
G(I)=7 
ENDIF 

IF((F(I).LT.E(I))AND.(F(I).LT.D(I))AND. 

+ (F(I).LT.P(I)))THEN 
G(I)=5 
ENDIF 
C ELSE 
C ENDIF 

IF((P(I) .LT.E(I)) AND.(P(I) .LT.D(I)) AND. 

+ (P(I).LT.F(I)))THEN 
G(I) = 3 
ENDIF 
C ELSE 
C ENDIF 

IF((D(I).EQ.E(I)).AND.(D(I).EQ.F(I)))THEN 

G(I)=9 

ENDIF 

IF((D(I).LT.F(I))AND.(E(I).LT.F(I)))THEN 

IF(D(I).EQ.E(I))THEN 

G(I)=4 

ENDIF 

ENDIF 

IF((D(I).LT.E(I))AND.(F(I).LT.E(I)))THEN 

IF(D(I).EQ.F(I))THEN 
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G(I)=6 
END IF 
ENDIF 

IF((F(I).LT.D(I))AND.(E(I).LT.D(I)))THEN 

IF(F(I).EQ.E(I))THEN 

G(I) = 8 

ENDIF 

ENDIF 

90 CONTINUE 

DO 1000 1 = 1,51 
DO 2000 J = 1,19 
PLOT(I,J) = G(J + 19*(I-1)) 

2000 CONTINUE 
1000 CONTINUE 

DO 3000 1 = 1^51 
C DO 4000 J = 1,42 

WRITE(8,5000)(PLOT(I,J),J = 1,19) 

3000 CONTINUE 
5000 format(20x,19il) 
stop 
end 
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1111111113111111337 

1111333111131153337 

1111111111111153537 

1135335115311111137 

1155151175711111117 

1111111177111777757 

1111111177111777777 

1155557777177777711 

1355555777711777717 

1111111111777777137 

1131111113377777117 

1131113111111111137 

1131351111131111737 

1131155357511117111 

1131355557711355311 

1111111333131133155 

1131131333131133155 

3113311353133315315 

3115111157111111117 

1355355157771137317 

1135555137771177117 

1131155111777777115 

1131135311177777135 

1135311315777777135 

1135313513577777135 

1131331353177177115 

1131331333137177711 

1111333771111117773 

1111155777111111175 

1355557777771177711 

1555557777771177735 

1531113777111777115 

1333133177711777111 

7315533577777771113 

7513551577777777775 

1351155317771171111 

1155555111777771113 

1135355337777777111 

1111151117771177111 

1355555777777777773 

1555555777777777775 

1151355777777777111 

1153553577777777711 


Best fit plot obtained for the cylinder belonging to set A. 
Numerals "1, 3, 5, T denote the original image, the 3 x 3 
image, the 5 x 5 image, and the 7x7 image respectively. 
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533333333333333 

333333333333333 

333333333333333 

333333333333333 

333333333333333 

333333333333333 

333333333333333 

333333333333333 

333333333333333 

333333333333333 

333333333333333 


Best fit plot for the sphere belonging to set B. 

Numerals "3, 5" denote the filtered 3x3 and 5x5 images 
of the original sphere. 
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55533331111111155555555 
55555555 1111 15555555555 
55555555555555555555555 
55555555555555555555555 
55555555555555555555555 
55555555555555555555555 
55555555555555555555555 
55555555555555555555555 
55555555555555555555555 
55555555555555555555555 
55555555555555555555555 
55555555555555555555555 
55555555555555555555555 
55555555555555555555555 
55555555555555555555555 
55555555555555555555555 
55555555555555555555555 
55555555555555555555555 
55555555555555555555555 
55555555555555555555555 
55555555555555555555555 


Best fit plot for the cylinder belonging to set B. 
Numerals "1, 3, 5" denote the original cylinder image, 
the 3 x 3 image, and the 5 x 5 image. 
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C***** PROGRAM SURFACE 

q ***************************************************** ********* 

C THIS PROGRAM APPROXIMATES THE COEFFICIENTS OF A SURFACE 
C GENERATED BY GIVEN DATA POINTS. THE INPUT FILE IS ’DATA.DAT’ 
C CONSISTING OF COORDINATES OF POINTS ON SOME SURFACE. 

C 

INTEGER U,K,IP 

REAL X(9000),Y (9000),Z(9000),X_2(9000) 

REAL Y 2(9000), Z_2(9000),P(9000,10) 

REAL YZ(9000),ZX(9000),XY (9000), P PTR(9000,10,10),SC(10,10) 

REAL A(4,4),B(6,4),B_TR(4,6),C(6,6),H(6,6),H_INV(6,6) 

REAL RIS(4,8)A_INV(4,4),BA_INV(6,4),BA_INVBT(6,6),M(6,6) 

REAL H_INVM(6,6),M PR(6,6),AI(6,6),BI(6,6),CI(6,6) 

REAL EIGVAL(6,6),EIGVEC(6,6),EI_VEC(6)A INVBT(4,6) 

REAL ALPHA(4),BETA(6)A_VECT(10) 

CHARACTER‘18 INFILE, OUTFILE 
TYPE*,’ ENTER COORDINATES FILE 
ACCEPT*, INFILE 

TYPE*,’ ENTER OUTPUT COEFFICIENTS FILE 
ACCEPT*, OUTFILE 

OPEN(UNIT = l.FILE = INF1LE, STATUS = ’OLD’) 

OPEN(UNIT = 2, FILE = OUTFILE, STATUS = ’NEW*) 

C* THE CONSTRAINT MATRIX H AND H INV IS CREATED ********* 

WRITE(*,3) 

3 FORMAT(5X,’INPUT TOTAL POINTS NOT EXCEEDING 7750: IP = ’) 
READ(*,*) IP 
ROOT = 1/(SQRT(2.)) 

DO 24 1 = 1,6 
DO 26 J = 1,6 
H(U) = 0 

26 CONTINUE 
24 CONTINUE 
H(l,l) = l 
H(2,2) = 1 
H(3,3) = 1 
H(4,4) = ROOT 
H(5,5) = ROOT 
H(6,6) = ROOT 
C 

ROOT1 = SQRT(2.) 

DO 20 1=1,6 
DO 22 J= 1,6 
H INV(U)=0 
22 CONTINUE 
20 CONTINUE 
H INV(1,1) = 1 
H'INV(2,2) = 1 
H"INV(33) = 1 
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H_INV(4,4) = ROOT1 
H INV(5,5) = ROOT1 
H~INV(6,6) = ROOT1 

C* DATA IS READ HERE **** ****** * 

DO 30 1 = 1, IP 
READ(1,*) (X(I),Y(I),Z(I)) 

30 CONTINUE 

C ****** THE VECTOR P FOR SCATTER MATRIX IS FORMED HERE ***** 

DO 32 1 = 1, IP 
X_2(I) = X(I)**2 
Y_2(I) = Y(I)**2 
Z 2(I) = Z(I)**2 
YZ(I) = Y(I)*Z(I) 

ZX(I) = Z(I)*X(I) 

XY(I)=X(I)*Y(I) 

32 CONTINUE 
DO 34 1 = 1, IP 
P(I,1) = X_2(I) 

P(I,2) = Y_2(I) 

P(U) = Z_2(I) 

P(I,4) = YZ(I) 

P(I,5) = ZX(I) 

P(I,6)=XY(I) 

P(I,7)=X(I) 

P(I,8) = Y(I) 

P(I,9) = Z(I) 

P(I,10) = 1 

34 CONTINUE 
DO 36 1 = 1, IP 
DO 38 J = 1,10 
DO 40 K=l,10 
P_PTR(I r I,K) = P(I r J)*P(I,K) 

40 CONTINUE 
38 CONTINUE 

36 CONTINUE 
DO 42 J = 1,10 
DO 44 K= 1,10 
SC(J,K)=0 
44 CONTINUE 

42 CONTINUE 

C**** THE SCATTER MATRIX IS FORMED HERE ************** 


DO 46 J= 1,10 
DO 48 K= 1,10 
DO 50 1 = 1, IP 

SC(J,K) = SC(J,K) + P PTR(U,K) 
50 CONTINUE 
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48 CONTINUE 
46 CONTINUE 

C******* THE SCATTER MATRIX SC IS DECOMPOSED INTO A,B,B_TR,C ** 

DO 52 I = 1,6 
DO 54 J = 1,6 
C(IJ)=SC(U) 

54 CONTINUE 

52 CONTINUE 
DO 56 1 = 1,6 
DO 58 J = l,4 
B(U)=SC(U+6) 

58 CONTINUE 

56 CONTINUE 
DO 60 1 = 1,4 
DO 62 J = l,6 
B_TR(I,J) = SC(I + 6,J) 

62 CONTINUE 

60 CONTINUE 
DO 64 1 = 1,4 
DO 66 J = l,4 
A(U)= SC(I+6,J+6) 

66 CONTINUE 

64 CONTINUE 
DO 68 1 = 1,4 
DO 70 J = l,4 
RIS(U)=A(U) 

70 CONTINUE 

68 CONTINUE 

CALL INVERS(RIS,4,4,8) 

DO 72 1 = 1,4 
DO 74 J = l,4 
A_INV(U) = RIS(U) 

74 CONTINUE 

72 CONTINUE 

C ***** NOW TO COMPUTE M 

DO 76 1=1,6 
DO 78 J = 1,4 
BA_INV(U)=0 
78 CONTINUE 

76 CONTINUE 
DO 80 1 = 1,6 
DO 82 J = 1,4 
DO 84 K= 1,4 

BA INV(I r I) = BA INV(I,J) + B(I,K)*A_INV(K r r) 

84 CONTINUE 

82 CONTINUE 

80 CONTINUE 
DO 86 1 = 1,6 
DO 88 J = 1,6 
BA_INVBT(U)=0 
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88 CONTINUE 
86 CONTINUE 
DO 90 1 = 1,6 
DO 92 J = 1,6 
DO 94 K= 1,4 

BA INVBT(U) = BA INVBT(U) + BA_INV(I,K)*B_TR(K r I) 
94 CONTINUE 

92 CONTINUE 

90 CONTINUE 
DO 96 1 = 1,6 
DO 98 J= 1,6 

M(U) = C(U)-BA_INVBT(U) 

98 CONTINUE 

96 CONTINUE 
C 

C ******** NOW TO COMPUTE M’ ******** * 

C 

DO 1001 = 1,6 
DO 102 J = 1,6 
H_INVM(U)=0 
102 CONTINUE 
100 CONTINUE 
DO 104 1=1,6 
DO 106 J = 1,6 
DO 108 K= 1,6 

H_INVM(U) = H_INVM(U) + H_INV(I,K)‘M(K r I) 

108 CONTINUE 
106 CONTINUE 
104 CONTINUE 
DO 110 1=1,6 
DO 112 J = 1,6 
M_PR(U)=0 
112 CONTINUE 
110 CONTINUE 
DO 114 1=1,6 


DO 116 J= 1,6 
DO 118 K= 1,6 

M PR(U) = M_PR(U) + H_INVM(I,K)*H_INV(K r I) 

118 CONTINUE 
116 CONTINUE 
114 CONTINUE 
C 

C **** NOW TO FIND THE EIGEN VALUES OF M’ ********** 

C 


ND=6 

CALL EIG(ND,M_PR,EIGVAL,EIGVEC) 


C 

C ******* TO FIND THE SMALLEST EIGEN VALUE AND ITS CORRESPONDING ** 
C ******* EIGENVECTOR * ********************* 

C 


S_EIG = EIGVAL(1,1) 
KOUNT = 1 
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DO 1201=2,6 

IF (S EIG.GT.EIGVAL(I,I)) THEN 
S_EIG = EIG VAL(I,I) 

KOUNT=I 

ENDIF 

120 CONTINUE 
DO 1221 = 1,6 

El VEC(I) = EIGVEC(I,KOUNT) 

122 CONTINUE 
DO 124 1 = 1,6 
BETA(I)=0 
DO 126 J = 1,6 

BETA(I) = BETA(I) + H_INV(U)*EI_VEC(J) 

126 CONTINUE 
124 CONTINUE 
DO 128 1 = 1,4 
DO 130 J = l,6 
A_INVBT(U)=0 
DO 132 K= 1,4 

A INVBT(U) = A_INVBT(U) + A INV(I,K)*B_TR(K r J) 

132 CONTINUE 
130 CONTINUE 
128 CONTINUE 
DO 1341 = 1,4 
ALPHA(I) = 0 
DO 136 J = 1,6 

ALPHA(I) = ALPHA(I) + A_INVBT(U)*BETA(J) 

136 CONTINUE 

ALPHA(I) = - ALPHA(I) 

134 CONTINUE 
DO 1381 = 1,6 
A_VECT(I) = BETA(I) 

138 CONTINUE 
DO 140 1 = 1,4 
A VECT(I + 6)=ALPHA(I) 

140 CONTINUE 
C DO 142 1=1,10 

WRITE(2,*) (’ THE INPUT FILE WAS INFILE, ’"’) 
WRITE(2,*) (’ THE OUTPUT FILE IS ”\OUTFILE,’"’) 
WRITE(2,*) (’ THE COEFF OF X-SQUARED IS ’,A_VECT(1)) 
WRITE(2,*) (’ THE COEFF OF Y-SQUARED IS ’/l_VECT(2)) 
WRITE(2,*) (’ THE COEFF OF Z-SQUARED IS ’A_VECT(3)) 
WRITE(2,*) (’ THE COEFF OF YZ IS ’,A_VECT(4)) 

WRITE(2,*) (’ THE COEFF OF ZX IS 'A VECT(5)) 

WRITE(2,*) (’ THE COEFF OF XY IS 'A~ VECT(6)) 

WRITE(2,*) (’ THE COEFF OF X IS >_VECT(7» 

WRITE(2,*) (’ THE COEFF OF Y IS ’,A_VECT(8)) 

WRITE(2,*) (’ THE COEFF OF Z IS 'A VECT(9)) 

WRITE(2,*) (’ THE CONSTANT D IS ’A_VECT(10» 

C142 CONTINUE 

CLOSE(UNIT = 2,DISPOSE = ’SAVE’) 

CLOSE(UNIT = l,DISPOSE = ’SAVE’) 
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END 


CC 

C * 


SUBROUTINE INVERS(RIS,N,NX,MX) 
DIMENSION RIS(NX,MX) 

N1 = N-1 
N2=2*N 
DO 2I = 1,N 
DO 1 J = 1,N 
J1=J + N 

1 RIS(I r Jl)=0. 

J1 = I + N 

2 RIS(U1) = 1. 

DO 10 K=1,N1 

C=RIS(K,K) 

IF (ABS(C)-O.OOOOOl) 3,3,5 

5 K1 = K+ 1 

DO 6 J = K1,N2 

6 RIS(K r I) = RIS(K r I)/C 
DO 10 I = K1,N 

C=RIS(I,K) 

DO 10 J = K1,N2 
RIS(U) = RIS(I r I)-C*RIS(KJ) 

10 CONTINUE 
NP1 = N+1 

IF (ABS(RIS(N,N))-0.000001) 3,3,19 

19 DO 20 J = NP1.N2 

20 RIS(N r I) = RIS(N r I)/RIS(N ) N) 

DO 200 L=1,N1 

K = N-L 
K1 = K+1 

DO 200 I=NP1,N2 
DO 200 J=K1,N 

200 RIS(K,I) = RIS(K,I)-RIS(KJ)*RIS(J,I) 

DO 250 I = 1,N 
DO 250 J = 1,N 
J1=J+N 

250 RIS(U) = RIS(IJ1) 

RETURN 

3 TYPE*, ’SINGULARITY IN ROW FOUND’ 
RETURN 

END 

SUBROUTINE EIG(ND/J,BI,CI) 

DIMENSION AI(ND,ND),BI(ND,ND),CI(ND,ND) 

INTEGER N1,M1,N2,M2 

N1 = ND 

M1 = ND 

N2=ND 

M2=ND 
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ANORM=0.0 

SN=FLOAT(N2) 

DO 100 I = 1,N2 
DO 101 J = 1,N2 
IF (I-J) 72,71,72 

71 BI(IJ) = 1.0 
GOTO 101 

72 BI(U)=0.0 

ANORM = ANORM + AI(U) *AI(U) 

101 CONTINUE 

100 CONTINUE 

ANORM = SORT (ANORM) 

FNORM = ANORM*(1.0E-09/SN) 

THR= ANORM 
23 THR=THR/SN 

3 IND=0 

DO 102 I = 2,N2 
11 = 1-1 

DO 103 J = 1,11 

IF (ABS(AI(J,I))-THR) 103,4,4 

4 IND = 1 
AL=-AI(J,I) 

AM = (AI(J,J)-AI(I,I))/2.0 

AO = AL/SQRT((AL*AL) + (AM* AM)) 

IF (AM) 5,6,6 

5 AO = - AO 

6 SINX = AO/SQRT(2.0*(1.0 + SQRT(1.0-AO*AO))) 
SINX2 = SINX*SINX 

COSX = SQRT(1.0-SINX2) 

COSX2=COSX*COSX 
DO 104 K=1,N2 
IF (K-J) 7,10,7 

7 IF (K-I) 8,10,8 

8 AT =AI(K r I) 

AI(K r J) = AT*COSX-AI(K,I)*SINX 
AI(K,I) = AT*SINX+AI(K,I)*COSX 
10 BT = BI(KJ) 

B^KJ) = BT*COSX-BI(K,I)*SINX 
BI(K,I) = BT*SINX + BI(K,I)*COSX 

104 CONTINUE 

XT = 2.0*AI(J,I)*SINX*COSX 
AT =AI(J r J) 

BT=AI(I,I) 

AI(J4) = AT*COSX2 + BT*SINX2-XT 

AI(I,I) = AT*SINX2 + BT*COSX2 + XT 

AI(J4) = (AT-BT)*SINX*COSX+AI(J,I)*(COSX2-SINX2) 

AI(U)=AI(J,I) 

DO 105 K= 1,N2 
AI(J,K)=AI(K r I) 

AI(I,K)=AI(K,I) 

105 CONTINUE 

103 CONTINUE 
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SS Kg 


102 CONTINUE 
IF (IND) 20,20,3 
IF (THR-FNORM) 25,25,23 
DO 110 I=2,N2 
J=I 

IF ((ABS(AI(J-l r J-l)))-(ABS(AI(J r I)))) 30,110,110 
AT=AI(J-l r I-l) 

AI(J-U-l)=AI(J r I) 

AI(J,J)=AT 
DO 111 K=1,N2 
AT=BI(K r I-l) 

BI(K r I-l) = BI(K r J) 

BI(ICJ)=AT 

111 CONTINUE 
J=J-1 

IF (J-l) 110,110,29 
110 CONTINUE 

DO 112 I = 1,N2 
DO 114 J = 1,N2 
ci(U)=Bi(y) 

BI(U) = AI(U) 

114 CONTINUE 

112 CONTINUE 
RETURN 

END 

C **************** ********** ******************************** 

c 
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